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ABSTRACT 

The angular power spectrum estimator developed by Peebles (1973) and 
Hauser & Peebles (1973) has been modified and applied to the 4 year maps 
produced by the COBE DMR. The power spectrum of the observed sky has 
been compared to the power spectra of a large number of simulated random 
skies produced with noise equal to the observed noise and primordial density 
fluctuation power spectra of power law form, with P(k) oc A;". The best fltting 
value of the spectral index in the range of spatial scales corresponding to spherical 
harmonic indices 3 < £ ^ 30 is an apparent spectral index Uapp = 1.131q;4 which 
is consistent with the Harrison-Zel'dovich primordial spectral index npri — I 
The best fitting amplitude for Uapp = I is (QjiMs)^'^ — 18 fjK. 
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1. Introduction 

The spatial power spectrum of primordial density perturbations, P{k) where k is the 
spatial wavenumber, provides evidence about processes occurring very early in the history 
of the Universe. In the first moments after the Big Bang, the horizon scale ct corresponds to 
a current scale that is much smaller than galaxies, so the assumption of a scale free form for 
P{k) on large scales is natural, which imphes a power law P{k) oc A;". The Poisson equation 

= AnGp implies that oc P{k)/k'^, and hence that the potential fluctuations 

on scale A oc 1/k are A^^ oc \(f){k)\'^k^ oc k'^'^. Harrison (1970), Zel'dovich (1972), and 
Peebles & Yu (1970) all pointed out that the absence of tiny black holes implies that 
the limit as ^ oo of Acj) must be small, so n ^ 1, while the large-scale homogeneity 
implied by the near isotropy of the Cosmic Microwave Background Radiation (CMBR) 
requires that the hmit as A; —> of A0 must be small, so n 1. Thus the prediction 
of a Harrison-Zel'dovich or n = 1 form for P(k) by an analysis that excludes all other 
possibilities is an old one. This particular scalc-frcc power law is scale-invariant because the 
perturbations in the metric (or gravitational potential) are independent of the scale. The 
inflationary scenario (Starobinsky 1980; Guth 1981) proposes a tremendous expansion 
of the Universe (by a factor ^ 10^^) during the inflationary epoch, which can convert 
quantum mechanical fluctuations on a microscopic scale during the inflationary epoch 
into Gpc-scale structure now. To the extent that conditions were relatively stable during 
the small part of the inflationary epoch which produced the Mpc to Gpc structures we 
now study, an almost scale-invariant spectrum is produced (Bardeen, Steinhardt & Turner 
1983). Bond & Efstathiou (1987) show that the expected variance of the coefficients a^m 
in a spherical harmonic expansion of the CMBR temperature given a power law power 
spectrum P{k) oc fc" is < a^^ > oc r[£ + (n - 1) /2]/r[i + (5 - n) /2] for i < 40. Thus a study 
of the angular power spectrum of the CMBR can be used to place limits on the spectral 
index n and test the inflationary prediction of a spectrum close to the Harrison-Zel'dovich 
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spectrum with n — 1. 

The angular power spectrum contains the same information as the angular correlation 
function, but in a form that simplifies the visualization of fits for the spectral index n. 
Furthermore, the off-diagonal elements of the covariance matrix have a smaller effect for 
the power spectrum than for the correlation function. However, with partial sky coverage 
the multipole estimates in the power spectrum are correlated, and this covariance must be 
considered when analyzing either the correlation function or the power spectrum. 

The power spectrum of a function mapped over the entire sphere can be derived easily 
from its expansion into spherical harmonics, but for a function known only over part of the 
sphere this procedure fails. Wright et al. (1994) have modified a power spectral estimator 
from Peebles (1973) and Hauser & Peebles (1973) that allows for partial coverage and 
applied this estimator to the DMR maps of CMBR anisotropy. We report here on the 
application of these statistics to the DMR maps based on all 4 years of data (Bennett et al. 
1996). Monte Carlo runs have been used to calculate the mean and covariance of the power 
spectrum. Fits to estimate {Q%ms)^'^ ^ maximizing the Gaussian approximation 
to the likelihood of the angular power spectrum are discussed in this paper. Since we 
only consider power law power spectrum fits in this paper, we use Q as a shorthand for 
(Qrms)^'^ or Qrms-ps, which is the RMS quadrupole averaged over the whole Universe, 
based on a power law fit to many multipoles. (QjiMs)'^'^ should not be confused with the 
actual quadrupole of the high galactic latitude part of the sky observed from the Sun's 
location within the Universe, which is the Qrms discussed by Bennett et al. (1992). 
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2. Estimating the Angular Power Spectrum 

Wright et al. (1994) have discussed the modification of the Hauser-Peebles angular 
power spectrum estimator for use on CMBR anisotropy maps. To allow for the cutting out 
of the galactic plane, new basis functions are defined using a modified inner product: 

</,>.?£^ (1) 

where j is an index over pixels, and Wj is the weight per pixel. In the galactic plane cut, 
Wj = 0. We have not used weights proportional to the number of observations, so wj = 1 
outside of the galactic plane cut. The custom galactic cut used in this paper basically 
follows sin|6| = 1/3 with extra cut added in Sco-Oph and Orion (Banday et al. 1996). A 
total of 3881 pixels are used, or 63% of the sky. 

The modified Hauser-Peebles method in Wright et al. (1994) used basis functions 
defined using 

^ J-, Fqq < FooFim > J-^ Fim' < Fim'-Pfm > /^.x 

f^im - i'tm 2^ ^ {1} 

< -f'OO-f'OO > m'=-l ^ ^Im'^lm' > 

where the are real spherical harmonics and the inner product < fg > is defined over the 
cut sphere. These functions G^rn ^-re orthogonal to monopole and dipole terms on the cut 
sphere. Call this the MD method since the basis functions are orthogonal to the monopole 
and dipole. Let the MDQ method use basis functions orthogonal to the monopole, dipole 
and quadrupole: 

< i*00-r00 > m'=-l ^ ^Im'-flm' > m'=-2 ^ ^2m'-f2m' > 

In this paper we have used the MDQ method so our results for ^ > 3 are completely 
independent of the quadrupole in the map. We have also tabulated the power in £ = 2 
which is computed with G2m'^ which are orthogonal to the monopole, dipole, and those 
components of the quadrupole which occur earlier in the sequence than m. Because the 
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galactic cut used is not a straight |6| cut, the different -F2m's are not quite orthogonal, and 
the definition of depends slightly on the ordering of the -F^m's- We use the ordering 
1, cos(f) cos 20, sin0, sin 20. With these basis functions we compute the power spectrum 
estimators 

2£ + 1 Yl,m=-i < ^'irrfi'im > 

which are quadratic functions of the maps. Note that for full sky coverage, is the 
variance of the sky in order £, but for partial sky coverage the response of Tj to inputs with 
V ^ i causes Tj to be larger than the order i sky variance. Table 1 of Wright et al. (1994) 
shows the input-output matrix for a straight 20° cut, while Table |l] shows the input-output 
matrix for the custom galaxy cut. The jump in Figure I at £ = 5 for the mean spectrum 
of Q = 17 /iK, n = \ inputs is caused by the off-diagonal response to £ = 3, while the 
off-diagonal response of £ = 4 to £ = 2 has been zeroed by the MDQ method. 

This method computes the power spectrum, a quadratic function of the map, which 
includes contributions from both the true sky signal and from instrument noise. We remove 
the contribution of the instrument noise by subtracting the power spectrum of a noise only 
map. This difference map can be constructed by subtracting the two maps made from the 
A and B sides of the DMR instruments: D = {A — B)/2. The sum map containing the 
real signal is S" = (A + B)/2. When we compute the quadratic power spectrum, we use 
the value 5*^ — = A x B, the power spectra reported here are the cross power spectra 
between the A and B sides of the DMR instrument. 

We have computed the power spectrum of the internal linear combination "free-free 
free" no galaxy (NG) map (Wright et al. 1994), Tng = -O.4512T31 + 1.2737r53 + O.3125T90, 
and the close to maximum signal-to-noise ratio maps based on O.6T53/O.93I -|- O.4T90/O.8I5. 
The denominators in the latter expression convert the Ray leigh- Jeans differential 
temperatures T53 and T90 into thermodynamic AT's, but this conversion is included in the 
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coefficients for T/vg- This process can also be applied to compute the cross power spectrum 
of the 53 GHz and 90 GHz maps by letting ^ = (53 + 90) and D = {53 - 90), after both 
maps have been converted into thermodynamic brightness temperature differences. Figure 
in shows the resulting power spectra for the three map combinations. 

We have binned the power spectra in quasi-logarithmic bins in i in order to minimize 
the increasing noise-to-signal ratio as i gets large. Figure ^ and Table ^ show these binned 
power spectral estimates. 



3. Monte Carlo Simulations 

In order to calibrate and test these methods for biases, it is necessary to simulate both 
the cosmic variance, which gives a random map with random spherical harmonic amplitudes 
chosen from a Gaussian distribution with a variance determined from the chosen Qj„ and 
riin, and the experimental variance, which gives the 360 million noise values needed per 
year. While programs to simulate the DMR time-ordered data do exist, we have not worked 
at this level of detail. Instead, we have used simulations that start with the maps. 

The effect of noise on the map production process can be simulated using 

T„ = a^A-'-'U (5) 

where ai is the noise in one observation, U is an uncorrelated vector of unit variance zero 
mean Gaussian random variables, and A is the matrix with diagonal elements An equal 
to the number of times the i^^ pixel was observed, and off-diagonal elements —Aij equal 
to the number of times the i^^ pixel was referenced to the j*^ pixel. Even though A is 
singular, Wright et al. (1994) give a rapidly convergent series technique for generating noise 
maps. Thus each noise map depends on 6144 independent Gaussian unit variance random 
variables and the parameter cxi. 
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The signal map that is added to the noise maps to give the "observed" maps is 
generated using independent Gaussian random amphtudes with variances given by Bond & 
Efstathiou (1987) for £ < 40. The simulations done here included €s up to 39, so the signal 
map depends on 1600 Gaussian independent unit variance random variables and the two 
parameters Qm and 

4. Mctximum Likelihood Estimation 

Using the Monte Carlo's, we find the mean power spectrum Tf{Qin,nin), and 
the covariance matrix C{Qin,nin)u' = {{Tf — TDiTfi — Tf)). For the actual power 
spectrum T/ from the real sky or a Monte Carlo simulation, define the deviation vector 
= Tf — T^iQimriin) and the statistic — e^C~^e. All of the fits in this paper are 
based on the range i = irain ■ ■ ■ ^max with imin = 3 and i^ax = 30. C is thus a 28 x 28 
matrix. Ignoring the quadrupole is reasonable because the galactic corrections are largest 
for ^ = 2, and the maximum order used is set by the DMR beam-size of 7° and the increased 
computer time required to analyze more orders. Since the magnitude of the covariance 
matrix gets larger rapidly when (5i„ increases there is a bias toward large values of Q when 
minimizing x^. One can allow for this by minimizing — 21n(L) instead of x^, where L is the 
Gaussian approximation to the likelihood: 

- 2 ln(L) = + ln(det(C)) + const. (6) 

For any given power spectrum, we can adjust Qm and n^^ until — 21nL is minimized. 
This gives us the maximum likelihood fit of a power law power spectrum to the given 
power spectrum. We have called the values of Q and n that maximize the likelihood for 
the observed power spectrum Qml and uml since these are maximum likelihood values. 
The maximum likelihood technique gives an asymptotically unbiased determination of the 
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amplitude Q and index n, but only as the observed solid angle goes to infinity. Since we 
are limited to about 8 sr of sky, asymptotically unbiased means biased in practice, both for 
the quadratic statistics considered here and for linear statistics used by Gorski et al. (1994) 
and Bond (1995). Our use of a Gaussian approximation to the likelihood for our quadratic 
statistics can introduce additional errors. We use our Monte Carlo simulations to calibrate 
our statistical methods to avoid biased final answers. 

Our k^^ Monte Carlo run depends on a set of random variables {Z^} (with 1600 + 
12288 elements for a cross-analysis needing two noise maps) having a known distribution, 
and the three parameters Qm-, rim and ai. ai can be determined with great precision using 
the time-ordered data. Hence one needs to run many Monte Carlo simulations with several 
different values for Qin and nj„ and compare the fitted values QML,k and nML,k to the 
fitted values for the real data, QML,ohs and nML,ohs- For the k^^ realization {Z^}, the fitted 
values QML,k and nML,k are a continuous function of the input parameters Qm and nj„, 
and one can choose values Qin = Qmatch,k and ni„ = rimatch^k such that QuL^k = QML,obs 
and nML,k = nML,obs- By choosing many different realizations of {Zk}, one creates many 
different {Q match, nmatch) pairs. Figure ^ shows this cloud of points for the 2 year 53 x 90 
cross-power spectrum, along with contours of — 21nL. The spectral index we give is 
the median of the set of rimatch^, and the 16%-tile to 84%-tile range in rimatch defines 
the ±1(7 range. The value l.lStoi given in the Abstract is the weighted mean of these 
determinations, but we have not reduced the error because the cosmic variance is common 
to all three maps. 

The value of (Q'rms)^'^ '^^^ found by doing a one parameter maximum likelihood 
fit for Qml with n fixed at 1. After using the Monte Carlo runs to debias the maximum 
likelihood results, we get the values shown in Table ^ We can also find the best fit values of 
Q for other values of n. The best fit Q values for n forced to be 1.25 are smaller than those 
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for n forced to be 1 by an amount that allows us to estimate the effective wavenumber of 
our amplitude determination. We find that i^ff = 7.9 for the 53+90 A x B case and 7.3 for 
the 53 X 90 case. In Figure |^ we have chosen to plot the amplitude at ^ = 7 which is closest 
to the effective wavenumber in order to minimize the correlation between the amplitude 
and the spectral index. 

5. Discussion 

The angular power spectrum of the four year COBE DMR maps has been calculated, 
and it is very consistent with a Harrison-Zel'dovich primordial spectrum Upri = 1, especially 
after the small correction for the "toe" of the Doppler peak which gives an expected 
apparent index of napp ~ 1.1 for Q = 1 CDM models. Models with a cosmological constant 
(ACDM) predict a smaller Uapp ~ 0.75 (Kofman & Starobinsky 1985) that is still consistent 
with the COBE DMR observations. The amplitude derived from this analysis is in between 
the {Q'rmsY'^ ~ 17 /iK derived from the first year maps and the {Q^msY'^ ~ ^9 A*-^ derived 
from the two year maps (Wright et al. 1994). The amplitude from the no galaxy map is 
consistent with but now slightly higher than the amplitude derived from the 53 and 90 
GHz maps, indicating that galactic contamination is not a major problem with the chosen 
galactic cut. 

We are grateful for the efforts of the COBE team and the support of the Office of 
Space Sciences at NASA. Charley Lineweaver provided helpful comments on an early draft 
of this paper. 
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Table 1. 10^ x the mean over m of the power spectra of F^/^ for the custom galaxy cut. 
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Table 2. Binned power spectra of the 4 yr DMR Maps. 



I Range 


4// 


NG >1 X S 


53 X 90 


53 + 90 >1 X S 


2 


2.1 


0.08 ±0.68 


0.16 ±0.65 


0.15 ±0.60 


3 


3.1 


0.76 ± 0.70 


0.99 ±0.59 


0.90 ±0.57 


4 


4.1 


1.65 ±0.76 


1.51 ±0.57 


1.52 ±0.56 


5-6 


4.6 


1.62 ±0.51 


1.40 ±0.35 


1.30 ±0.34 


7-9 


6.3 


1.71 ±0.52 


0.85 ± 0.29 


0.94 ±0.28 


10-13 


8.9 


-0.02 ± 0.79 


0.78 ± 0.33 


1.19 ±0.30 


14-19 


11.5 


1.50 ± 1.42 


1.75 ±0.46 


1.52 ±0.40 


20-30 


13.2 


0.37 ±4.81 


0.57 ± 1.39 


0.02 ±1.08 



Note. — All values have been normalized to the mean for 
Q = 17 jjK, n — 1 Monte Carlo runs. 



Table 3. Power Law Fits to 4 year DMR map power spectra. 



Maps Happ {QuMsf''' at napp = 1, [/iK] 



53+90 AxB l.nt^fi 17.91}:^ 

53x90 l-0'2toA7 17.3t2:? 
NG AxB 1.22l[J:?? 19.6lij 
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Fig. 1.— Cross power spectra for the 53+90 A x 53 x 90 and NG ^ x S maps. 
measures the variance of the sky due to order £ harmonics for full sky coverage, but partial 
sky coverage changes the expected value slightly as seen in the curves showing the average 
power spectrum oi Q — 17 /iK, n — 1 models in the cut sky. Values are shifted upward by 
400 for NG and 900 for 53 x 90, as shown by the horizontal lines marking zero power. 

Fig. 2. — Binned cross power spectra for the 53+90 A x B, 53 x 90 and NG Ax B maps, 
normalized to the mean power spectrum of Q = 17 fiK, n = 1 simulations, plotted on a 
logarithmic scale, igff is the effective wavenumber of the bin for n = 1. The thin curves 
show a CDM model with Upri — 0.96 including the effect of gravitational waves derived from 
Crittenden et al. (1993). 

Fig. 3. — Each point is an input parameter set that is consistent with the real 4 year 53 x 90 
DMR data for a given realization of the random cosmic and radiometer variance processes. 
The amplitude is specified using the RMS AT due to £ = 7 spherical harmonics because i^ff 
for this fit is 7.3. The hkelihood contours are at A(— 21nL) = 1, 4 and 9. 
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